Deterministic models
Forrest W. Crawford
Overview
We begin with deterministic compartmental models, which are usually straightforward to specify. With the righ software, they are straightforward to solve. These models find application in virtually every scientific field. In public health, we see them most often in infectious disease epidemiology, cancer biology, pharmacokinetics, and healtcare operations research. We will use examples from Brauer, Castillo-Chavez, and Castillo-Chavez (2001 and others).
also here: https://institutefordiseasemodeling.github.io/Documentation/general/model-seir.html
The focus here is on specifying systems of ordinary differential equations, not on solving those systems analytically. There is a huge literature on analytic solution – you have probably seen some of this in a calculus or differential equations course – but it’s beyond the scope of this course, which is focused on the modeling part. We will also not focus on qualitative dynamics and equilibrium analysis of DE systems. Instead, we will focus on specification and construction of systems, and using a computer to solve them numerically.
Examples of deterministic models
The most boring model (part 1)
Consider the most boring model for a quantity \(y(x)\) \[ \frac{dy}{dx} = \beta \] with initial condition \(y(0)=\alpha\). The rate of change of \(y\) as a function of \(x\) is constant.
What function of \(t\) has slope \(\beta\) for every \(t\)? The linear function \[ y(x) = \alpha + \beta x . \] Have you seen this before?
The most boring model (part 2)
Earlier we constructed the linear “regression” model \[ y = \alpha + \beta x \] from a differential equation.
If \(x\) is a treatment, and we regard \(y(x)\) as the potential outcome when treatment is set to \(x\), then the “treatment effect” is the change in \(y\) induced by a unit change in \(x\), \[ \frac{dy}{dx} = \beta \] This is one reason that people (economists mostly) think of causal effects as derivatives of outcome functionals.
Birth process
Recall the basic birth process, where the rate of growth in \(y(t)\) is proportional to \(y(t)\) itself, \[ \frac{dy}{dt} = \lambda y(t) \] with initial condition \(y(0)=y_0\), which has solution \[ y(t) = y_0 e^{\lambda t} . \] When \(\lambda>0\), \(y(t)\) increases “exponentially”.
Death process
The “death” process is the same, but we let \(\lambda=-\mu<0\). Then \[ \frac{dy}{dt} = -\mu y(t) \] with initial condition \(y(0)=y_0\), which has solution \[ y(t) = y_0 e^{-\mu t} \] and \(y(t)\) decreases/decays “exponentially” toward zero.
Logistic growth I
Recall the exponential growth model \[ \frac{dy}{dt} = \lambda y(t) \] with initial condition \(y(0)=y_0\). Its solution is \(y(t) = y_0 e^{\lambda t}\). As \(t\to\infty\), \(y(t)\to\infty\). This unbounded growth may be unrealistic. For example, resource constraints may prevent cells from dividing indefinitely.
Suppose \(K\) is the “carrying capacity” of the medium, and \(\lambda\) is the intrinsic gowth rate as above. Consider the model \[ \frac{dy}{dt} = \lambda y(t)\left( 1 - \frac{y(t)}{K}\right) \] with initial condition \(y(0) = y_0\). This is just like the pure birth model, except there is a negative quadratic term.
Logistic growth II
\[ \frac{dy}{dt} = \lambda y(t)\left( 1 - \frac{y(t)}{K}\right) \]
\(y(t)\) is constant when \(dy/dt=0\), which occurs whe either \(y(t)=0\) or \(y(t)=K\). Therefore
- If \(y(t) < K\), then \(dy/dt>0\) and the population size increases toward \(K\).
- If \(y(t) > K\) then \(dy/dt<0\) and the population size decreases toward \(K\).
In fact, the solution can be written down precisely: \[ y(t) = \frac{K y_0}{y_0 + (K-y_0)e^{-\lambda t}} \]
ts = seq(0,10,by=0.01)
lam = 0.5
K = 0.8
y0s = c(0.05, 0.2,0.5,1)
plot(0,
type="n",
xlim=range(ts),
ylim=c(0,1),
ylab="Number of bacteria y(t)",
xlab="Time t",
bty="n")
text(10, K, "K=0.8", pos=3)
abline(h=K, lty="dashed", col="gray")
for(y0 in y0s) {
ys = K * y0 / (y0 + (K-y0)*exp(-lam*ts))
lines(ts, ys)
text(1,ys[100], bquote(y[0] == .(y0)), pos=1)
}
Logistic growth: SI model
Consider a population of individuals where every individual is susceptible or infected with an infectious disease. Let \(S(t)\) be the proportion of susceptibles, and let \(I(t)\) be the proportion of infectives. Then the system is characterized by \[ \frac{dI}{dt} = \beta S(t) I(t) \] where \(S(t) + I(t) = 1\). We can rewrite this model as a logistic growth model \[ \frac{dI}{dt} = \beta I(t)(1-I(t)) \] where initially we have \(I(0)=i_0 \le 1\) infectives. The carrying capacity is \(K=1\) and \(\lambda=\beta\) and \(y_0=I(0)\) above. The difference is that \(K=1\) is the maximum population proportion, so it is not possible to start at values of \(i_0\) greater than 1.
SIR model
Consider a population of \(N\) individuals where every individual is susceptible or infected or recovered. Let \(S(t)\) be the number of susceptibles, \(I(t)\) be the number of infectives, \(R(t)\) the number of recovered. Then the system is characterized by
\[ \begin{aligned} \frac{dS}{dt} &= -\beta S(t) I(t) \\ \frac{dI}{dt} &= \beta S(t) I(t) - \gamma I(t) \\ \frac{dR}{dt} &= \gamma I(t) \end{aligned} \]
with initial condition \(S_0\), \(I(0)\), \(R(0)\) and conservation equation \(S(t)+I(t)+R(t)=N\). This system does not have a solution that is easy to write down, but it is easy to understand!
SIR model with seasonal forcing
Suppose the transmission rate is a function of time, \(\beta(t) = a(1 + \sin(b + t/c))\)
\[ \begin{aligned} \frac{dS}{dt} &= -\beta(t) S(t) I(t) \\ \frac{dI}{dt} &= \beta(t) S(t) I(t) - \gamma I(t) \\ \frac{dR}{dt} &= \gamma I(t) \end{aligned} \]
When might periodicity, or seasonal forcing, be appropriate?
SEIR model
Now suppose there is a latent period following infection that delays onset of infectiousness. Individuals in this latent state are exposed, but not yet infectious.
\[ \begin{aligned} \frac{dS}{dt} &= -\beta S(t) I(t) \\ \frac{dE}{dt} &= \beta S(t) I(t) - \delta E(t) \\ \frac{dI}{dt} &= \delta E(t) - \gamma I(t) \\ \frac{dR}{dt} &= \gamma I(t) \end{aligned} \]
For what infectious diseases might a latent period be important? How does this change our interpretation of \(I(t)\)?
SEIR model with vital dynamics
Let’s incorporate births and deaths. Now the population is no longer closed, so let’s model the absolute number of individuals in each compartment instead of the proportion. Let \(N(t) = S(t) + E(t) + I(t) + R(t)\). Let \(\nu\) be the per-person birth rate, and let \(\nu\) be the per-person death rate.
\[ \begin{aligned} \frac{dS}{dt} &= \mu N(t) - \nu S(t) -\beta S(t) I(t) \\ \frac{dE}{dt} &= \beta S(t) I(t) - (\nu + \delta) E(t) \\ \frac{dI}{dt} &= \delta E(t) - (\nu + )\gamma I(t) \\ \frac{dR}{dt} &= \gamma I(t) - \nu R(t) \end{aligned} \]
The biggest difference in this model is that an “outbreak” can re-start after slowing down due to re-introduction of new susceptibles. When are vital dynamics important?
SEIRS (loss of immunity)
\[ \begin{aligned} \frac{dS}{dt} &= -\beta S(t) I(t) + \sigma R(t) \\ \frac{dE}{dt} &= \beta S(t) I(t) - \delta E(t) \\ \frac{dI}{dt} &= \delta E(t) - \gamma I(t) \\ \frac{dR}{dt} &= \gamma I(t) - \sigma R(t) \end{aligned} \]
When is loss of immunity due to prior infection important? For which diseases?
Lotka-Volterra predator-prey dynamics (optional)
Let \(x(t)\) be the number of prey, and let \(y(t)\) be the number of predators.
- Having abundant food, prey reproduce naturally in proportion (\(\alpha\)) to their number
- Prey die (are eaten) in proportion (\(\beta)\) to the product of their abundance and the number of predators
- Predators reproduce in proportion (\(\delta\)) to the product of their number and the number of prey they can eat
- Predators die in proportion to their number
\[ \frac{dx}{dt} = \alpha x(t) - \beta x(t) y(t) \] \[ \frac{dy}{dt} = \delta x(t) y(t) - \gamma y(t) \] with initial conditions \(x(0)\) and \(y(0)\).
Example: HIV care pipeline
\[ \begin{aligned} \frac{dD}{dt} &= \lambda_D I(t) - (\lambda_L + \mu) D(t) \\ \frac{dL}{dt} &= \lambda_L D(t) - (\lambda_A + \mu) L(t) \\ \frac{dA}{dt} &= \lambda_A L(t) - (\lambda_S + \mu) A(t) \\ \frac{dS}{dt} &= \lambda_S A(t) - \mu S(t) \\ \end{aligned} \]
Adapted from Gonsalves et al. (2017).
Example: HIV care pipeline
Gonsalves et al. (2017)
Drug use initiation and treatment
Let \(S(t)\) be the proportion of the population who do not use opiates, \(D(t)\) the population that do use opiates, and \(T(t)\) the proportion in treatment.
\[ \begin{aligned} \frac{dS}{dt} &= -\alpha S(t) -\beta S(t) D(t) + \sigma R(t) \\ \frac{dD}{dt} &= \delta E(t) - \gamma D(t) \\ \frac{dT}{dt} &= \gamma D(t) - \sigma T(t) \end{aligned} \]
cite olya’s and margret’s and ed’s paper?
Pharmacokinetic models
Let \(G(t)\) be the amount of drug in the GI tract, \(B(t)\) in the bloodstream, and \(O(t)\) in the organ of interest. Let \(I(t)\) be the known input dynamics of drug administration (oral, IV, etc).
\[ \begin{aligned} \frac{dG}{dt} &= I(t) - \lambda_G G(t) \\ \frac{dB}{dt} &= \lambda_G G(t) - \mu_B(t) - \lambda_T B(t) \frac{dO}{dt} &= \lambda_T B(t) - \mu_O O(t) \end{aligned} \]
Solving ODEs numerically
ODE solvers
Lots of ODEs can’t be solved analytically. For (well-behaved) systems that have a unique solution, it is often possible to solve for \(y(t)\) numerically, without having any calculus-style analytic insight into the structure of the solution. There are very fast computer algorithms for doing this.
We will focus on methods implemented in standard software packages.
General-purpose ODE solver: lsoda
library(deSolve)you will need to install this.
For help, look at ?deSolve
Example: SIR solution
library(deSolve)
sir.model <- function (t, x, params) {
S <- x[1]
I <- x[2]
R <- x[3]
with( as.list(params), {
dS <- -beta*S*I
dI <- beta*S*I-gamma*I
dR <- gamma*I
dx <- c(dS,dI,dR)
list(dx)
})}
times <- seq(0,120,by=1)
params <- c(beta=0.4,gamma=1/7)
xstart <- c(S=9999/10000,I=1/10000,R=0)
out <- as.data.frame( lsoda(xstart, times, sir.model, params) )
plot(out$time, out$I, ylab="Population proportion", xlab="Time",
type='l', bty="n", ylim=c(0,1), col="red")
lines(out$time, out$S, col="green")
lines(out$time, out$R, col="blue")
legend(100,0.5, c("S(t)", "I(t)", "R(t)"), lty=1, col=c("green", "red", "blue"))
Example: Lotka-Volterra solution
library(deSolve)
lv.model <- function (t, x, params) {
X <- x[1]
Y <- x[2]
with( as.list(params), {
dX <- alpha*X - beta*X*Y
dY <- delta*X*Y - gamma*Y
dx = c(dX,dY)
list(dx)
})}
times <- seq(0,120,by=1)
params <- c(alpha=1.1,beta=0.4,delta=0.1,gamma=0.4)
xstart <- c(X=0.5,Y=0.6)
out <- as.data.frame( lsoda(xstart, times, lv.model, params) )
plot(out$time, out$X, ylab="Population proportion", xlab="Time",
type='l', bty="n", col="red")
lines(out$time, out$Y, col="blue")
legend(100,0.9*max(out$X), c("X(t)", "Y(t)"), lty=1, col=c("red", "blue"))
Step 1: specify the model
sir.model <- function (t, x, params) {
S <- x[1]
I <- x[2]
R <- x[3]
with(as.list(params), {
dS <- -beta*S*I
dI <- beta*S*I-gamma*I
dR <- gamma*I
dx <- c(dS,dI,dR)
list(dx)
})}First, we see that we have created a function sir.model, which is a funciton of time t, a model state x, and parameters params. The state x is a three-component vector corresponding to the S, I, and R states of the model.
Within the function definition, we decompose x into S, I, and R states and store these separately for later use.
The with function takes a list of parameters and does an action (contained in the curly brace block) using those parameters in memory. The action here is definition of the derivatives dS, dI, and dR, and putting these together in a list.
Keep in mind that R is not doing anything fancy here. No calculus. At this point your computer does not know what an ODE system or SIR model is. We have just created a function that takes some arguments and returns threee numbers.
Set up the times, parameters, and starting values
times <- seq(0,120,by=1)
params <- c(beta=0.4,gamma=1/7)
xstart <- c(S=9999/10000,I=1/10000,R=0)Running lsoda
Here is the most important part
out <- as.data.frame( lsoda(xstart, times, sir.model, params) )The lsoda function takes these arguments and solves the system for every time in times. We store the result as a data frame instead of as a list.
Here is what out looks like
head(out)## time X Y
## 1 0 0.500000 0.6000000
## 2 1 1.224622 0.4358909
## 3 2 3.148179 0.3580257
## 4 3 8.168257 0.4066573
## 5 4 19.243069 1.0178402
## 6 5 17.947209 6.2733333Overview: how to model with ODEs
ODEs made easy
- Write the components (“compartments” or “types” or “states”) in the time-dependent system you want to model.
- Write the time-derivatives of the components as a function of time, and other compartments.
- Write the initial conditions for all the components
- Simulate the dynamics of the system
ODEs are most interesting when there are interactions between compartments.
A brief caveat
Modeling with ODEs seems easy – we just write down the rates of the system, then solve it numerically. For most easy-to-specify models you will encounter in public health research, this is a good plan. But it is sometimes possible to specify an ODE system whose solution does not exist.
Consider a system
\[ \frac{d y}{dt} = F(y,t) . \]
If \(F(y,t)\) is not continuous, a solution may not exist, or it may not be unique.
References
Brauer, Fred, Carlos Castillo-Chavez, and Carlos Castillo-Chavez. 2001. Mathematical Models in Population Biology and Epidemiology. Vol. 40. Springer.
Gonsalves, Gregg S, A David Paltiel, Paul D Cleary, Michael J Gill, Mari M Kitahata, Peter F Rebeiro, Michael J Silverberg, et al. 2017. “A Flow-Based Model of the Hiv Care Continuum in the United States.” Journal of Acquired Immune Deficiency Syndromes (1999) 75 (5). NIH Public Access: 548–53.